################################################
# Analysis for:
#
# Haglin, Kathryn, Soren Jordan, and Grant Ferguson. ``They’re Coming For You! 
#  How Perceptions of Automation Affect Public Support for Universal Basic Income''
#  Social Science Computer Review.
#  DOI: 10.1177/08944393231212252
#
# Required files: UBI automation.csv
#
################################################

library(tidyverse)

data <- read_csv("UBI automation.csv")
names(data)

# IDEOLOGY:
#	 1 = Extremely/very conservative
#	 2 = Conservative
#	 3 = Slightly conservative
#	 4 = Moderate, middle of the road
#	 5 = Slightly liberal
#	 6 = Liberal
#	 7 = Extremely/very liberal

# PID:
#	1 = Strong Republican
#	2 = Not very strong Republican/Don't know how strong (but Republican)
#	3 = Lean Republican
#	4 = Neither/true independent/don't kow
#	5 = Lean Democrat
#	6 = Not very strong Democrat/Don't know how strong (but Democrat)
#	7 = Strong Democrat

# INDUSTRY_CODE: industry matched to BLS code by authors

# All variables ending in _AUTO: how likely is it that the respective job will be automated?
#	1 = Not at all likely
#	2 = Not very likely
#	3 = Somewhat likely
#	4 = Very likely
#	5 = Don't know

# ECON_VULNERABLE: Economic vulnerability on a 0-100 scale (slider)

# UBI_SUPPORT:
#	1 = not at all likely to support
#	2 = not very likely to support
#	3 = neither likely nor unlikely
#	4 = somewhat likely to support
#	5 = very likely to support

# DEFICIT:
# 	1 = Extremely important (to reduce the deficit)
# 	2 = Very important
# 	3 = Moderately important
# 	4 = Slightly important
# 	5 = Not at all important (to reduce the deficit)

# INCOME:
# 	1 = Less than $10,000
#	2 = $10,000 - $19,999
#	3 = $20,000 - $29,999
#	4 = $30,000 - $39,999
#	5 = $40,000 - $49,999
#	6 = $50,000 - $59,999
#	7 = $60,000 - $69,999 
#	8 = $70,000 - $79,999
#	9 = $80,000 - $89,999
#	10 = $90,000 - $99,999
#	11 = $100,000 - $149,999
#	12 = More than $150,000

# EDUC:
# 	1 = Less than 6th grade
#	2 = Did not graduate high school
#	3 = High school equivalent
#	4 = Some college but no degree
#	5 = Associate's degree
#	6 = Bachelor's degree
#	7 = Advanced degree

# TREATMENT (unrelated survey experiment: DOI: 10.1111/SPOL.12760):
#              VAL
#      N        C       P
#    _______________________
# N |   1   |   2   |   3   |		
#   |_______|_______|_______|
# C |   4   |   0   |   6   |		POL
#   |_______|_______|_______|
# P |   5   |   7   |   8   |		
#   |_______|_______|_______|

# FEMALE:
#	1 = Female
#	0 = Non-female 

# RETRO:
# 	1 = Gotten much better (in the last 12 months, the economy ....)
# 	2 = Gotten somewhat better
# 	3 = Stayed about the same 
# 	4 = Gotten somewhat worse
# 	5 = Gotten much worse

# PROSPECT:
# 	1 = Get much better (in the last 12 months, the economy ....)
# 	2 = Get somewhat better
# 	3 = Stay about the same 
# 	4 = Get somewhat worse
# 	5 = Get much worse

# VOTEYES:
#	1 = Voted for president in 2016
#	0 = Non-female 

# Create binary versions of ideology and party ID

data$IDEO_BIN <- data$IDEOLOGY - 4
data$IDEO_BIN[data$IDEO_BIN < 0] <- -1
data$IDEO_BIN[data$IDEO_BIN > 0] <- 1
table(data$IDEOLOGY, data$IDEO_BIN)

data$PID_BIN <- data$PID - 4
data$PID_BIN[data$PID_BIN < 0] <- -1
data$PID_BIN[data$PID_BIN > 0] <- 1
table(data$PID, data$PID_BIN)

# Remove NA and SAHM/RETURED/Student/vague from industry code: we are 
#  interested in workers in the workforce and their perceptions
data$INDUSTRY_CODE_SLIM <- data$INDUSTRY_CODE
data$INDUSTRY_CODE_SLIM[data$INDUSTRY_CODE_SLIM == 0] <- NA
data$INDUSTRY_CODE_SLIM[data$INDUSTRY_CODE_SLIM > 70] <- NA

# Baseline: self-employed
data$INDUSTRY_CODE_SLIM_BASE70 <- data$INDUSTRY_CODE_SLIM
data$INDUSTRY_CODE_SLIM_BASE70[data$INDUSTRY_CODE_SLIM_BASE70 == 70] <- 1


length(na.omit(data$INDUSTRY_CODE_SLIM_BASE70))


# Make mini-dataset of industry (BLS code), industry (name), McKinsey score (see the
#  McKinsey study in the Appendix), and Low/Mid/High score. If there is no McKinsey score,
#  just record the industry (from BLS)

industry.auto.scores <- data.frame(code = 11, Industry = "Management", McKinsey = 35, Threat = "Low", numthreat = 1)
	
industry.auto.scores <- rbind(industry.auto.scores, 
			c(13, "Business/Finance", 43, "Medium", 2),
			c(15, "Computer/Math", 36, "Low",  1),
			c(17, "Architecture/Engineering", NA, NA, NA),
			c(19, "Life/Physical/Social Sciences", NA, NA, NA),
			c(21, "Community Services", 36, "Low", 1),
			c(23, "Legal", NA, NA, NA),
			c(25 ,"Education", 27, "Low", 1),
			c(27, "Arts/Entertainment", 41, "Medium", 2),	
			c(29, "Healthcare", 36, "Low", 1),
			c(31, "Healthcare Support", 36, "Low", 1),
			c(33, "Protective", 49, "Medium", 2),	
			c(35, "Food Prep", 73, "High", 3),
			c(37, "Building/Maintenance", NA, NA, NA),
			c(39, "Personal Care/Services", NA, NA, NA),
			c(41, "Sales", 53, "High", 3),	
			c(43, "Administrative", 39, "Low", 1),	
			c(45, "Farming/Fishing", 57, "High",3),	
			c(47, "Construction", 47, "Medium", 2),
			c(49, "Installation/Repair", 44, "Medium", 2),	
			c(51, "Production", 60, "High", 3),	
			c(53, "Transportation", 60, "High", 3),	
			c(55, "Military", 35, "Low", 1),
			c(60, "Freelance", NA, NA, 4),
			c(70, "Self-employed", NA, NA, 5))


# Add these low/medium/high scores to the survey dataset (to the industries which have a close match)
data$YOUR_JOB_AUTO_LOW_MID_HIGH <- NA
for(i in 1:length(data$INDUSTRY_CODE_SLIM)) {
	data$YOUR_JOB_AUTO_LOW_MID_HIGH[i] <- ifelse(data$INDUSTRY_CODE_SLIM[i] %in% as.numeric(industry.auto.scores$code), 
												as.numeric(industry.auto.scores$numthreat)[data$INDUSTRY_CODE_SLIM[i] == industry.auto.scores$code], 
												data$YOUR_JOB_AUTO_LOW_MID_HIGH[i])
}

table(data$YOUR_JOB_AUTO_LOW_MID_HIGH, data$INDUSTRY_CODE_SLIM)





####################################################
# Figure 1. Where do people work?

perc.v.real.part <- data %>%
	group_by(INDUSTRY_CODE_SLIM, PID_BIN) %>%
	count() %>%
	pivot_wider(names_from = PID_BIN, values_from = n)

full.count <- data %>%
	group_by(INDUSTRY_CODE_SLIM) %>%
	count()

perc.v.real.part$Industry <- NA

for(i in 1:length(perc.v.real.part$Industry)) {
	perc.v.real.part$Industry[i] <- ifelse(perc.v.real.part$INDUSTRY_CODE_SLIM[i] %in% industry.auto.scores$code,
									industry.auto.scores$Industry[industry.auto.scores$code == perc.v.real.part$INDUSTRY_CODE_SLIM[i]],
									perc.v.real.part$Industry[i])
}

# Rename variables, lengthen for ggplot
perc.v.real.part$Republicans <- perc.v.real.part$"-1"
perc.v.real.part$Independents <- perc.v.real.part$"0"
perc.v.real.part$Democrats <- perc.v.real.part$"1"
perc.v.real.part$Full <- full.count$n

perc.v.real.part <- perc.v.real.part %>%
	pivot_longer(cols = Republicans:Full,
				names_to = "Sample")

# Select only industry named variables, remove NA for plotting
perc.v.real.part <- perc.v.real.part %>%
	select(Industry, Sample, value) %>%
	na.omit()

# What's the numeric cutoff for the top 10 industries in the full sample?
top10.cutoff <- as.numeric(names(table(perc.v.real.part$value[perc.v.real.part$Sample == "Full"])[length(names(table(perc.v.real.part$value[perc.v.real.part$Sample == "Full"]))) - 10]))

top10 <- perc.v.real.part$Industry[perc.v.real.part$Sample == "Full" & perc.v.real.part$value > top10.cutoff]

# Make small dataset of just the top ten industries
perc.v.real.part.small <- perc.v.real.part[perc.v.real.part$Industry %in% top10 & !(perc.v.real.part$Sample %in% "Independents"),]

# Figure 1. Industry by partisanship.
ggplot(data = perc.v.real.part.small, aes(x = Sample, y = value, fill = Industry)) +
	geom_bar(stat = "identity", position = "dodge", color = "black") +
	theme_bw() +
	# PuRd color scale with 10th color of almost black added
	scale_fill_manual(values = c("#F7F4F9", "#E7E1EF", "#D4B9DA", "#C994C7", "#DF65B0", 
		"#E7298A", "#CE1256", "#980043", "#67001F", "#1a050d")) +
	ylab("N")
dev.off()

####################################################




####################################################
# Figure 2. How do perceptions compare with reality?
# Lack of relationship between actual threat and perception
# Democrat and Republican perceptions of threat
#	show that this difference doesn't really exist
# Does this track with economic vulnerability?

perc.v.real <- data %>% 
	group_by(INDUSTRY_CODE_SLIM) %>%
	summarise(Fear = mean(YOUR_JOB_AUTO, na.rm = TRUE),
				N = length(YOUR_JOB_AUTO))


perc.v.real$Industry <- NA
perc.v.real$McKinsey <- NA
perc.v.real$Threat <- NA

for(i in 1:length(perc.v.real$Industry)) {
	perc.v.real$Industry[i] <- ifelse(perc.v.real$INDUSTRY_CODE_SLIM[i] %in% industry.auto.scores$code,
									industry.auto.scores$Industry[industry.auto.scores$code == perc.v.real$INDUSTRY_CODE_SLIM[i]],
									perc.v.real$Industry[i])
	perc.v.real$McKinsey[i] <- ifelse(perc.v.real$INDUSTRY_CODE_SLIM[i] %in% industry.auto.scores$code,
									industry.auto.scores$McKinsey[industry.auto.scores$code == perc.v.real$INDUSTRY_CODE_SLIM[i]],
									perc.v.real$McKinsey[i])
	perc.v.real$Threat[i] <- ifelse(perc.v.real$INDUSTRY_CODE_SLIM[i] %in% industry.auto.scores$code,
									industry.auto.scores$Threat[industry.auto.scores$code == perc.v.real$INDUSTRY_CODE_SLIM[i]],
									perc.v.real$Threat[i])
								
}

# Remove healthcare support as it overplots
perc.v.real$Industry[perc.v.real$Industry == "Healthcare Support"] <- NA

perc.v.real <- na.omit(perc.v.real)

perc.v.real$Threat <- factor(perc.v.real$Threat, levels = c("Low", "Medium", "High"))
perc.v.real$McKinsey <- as.numeric(perc.v.real$McKinsey)




# Figure 2. fear of automation versus actual likelihood
ggplot(data = perc.v.real, aes(x = Fear, y = McKinsey, size = N, label = Industry, col = Threat)) +
	geom_label() + 
	scale_color_manual(values = c("#91003F", "#E7298A", "#C994C7")) +
	xlim(1.25, 2.75) + 
	guides(scale = "none") + 
	theme_bw() +
	theme(legend.position = c(0.2, 0.80)) +
	labs(x = "Average Self-Assessed Likelihood of Own Job Being Automated (by Industry)\n(1 = Not at all likely, 4 = Very likely)", 
		y = "Actual Threat of Automation (Expert Evaluation)", 
		col = "Actual Threat of Automation\n(Three Category)", size = "N in Industry in Full Sample") 
dev.off()

####################################################





####################################################
# Figure 3. What industries are perceived to be vulnerable? (full sample)

data.barplots.int <- data %>% 
				select(ends_with("AUTO")) %>%
				pivot_longer(cols = ends_with("AUTO"), names_to = "Occupation")

data.barplots.int$value[is.na(data.barplots.int$value)] <- 5
		
data.barplots <- data.barplots.int %>%
				mutate(Occupation = fct_relevel(fct_recode(factor(Occupation), 
					"Fast Food" = "FAST_FOOD_AUTO", "Insurance Claims" = "INSURE_CLAIM_PROCESS_AUTO", 
					"Software Engineer" = "SOFTWARE_ENGINEER_AUTO", "Legal Clerk" = "LEGAL_CLERK_AUTO",
					"Construction Work" = "CONS_WORK_AUTO", "Teacher" = "TEACHER_AUTO",
					"Nurse" = "NURSE_AUTO", "Politician" = "POLITICIAN_AUTO",
					"Your Job" = "YOUR_JOB_AUTO"))) %>%
				mutate(Rating = fct_relevel(fct_recode(factor(value), 
					"Not at all likely" = "1",
					"Not very likely" = "2",
					"Somewhat likely" = "3",
					"Very likely" = "4",
					"Don't know" = "5")))				

ggplot(data.barplots, aes(x = Occupation, fill = Rating)) +
	geom_bar(stat = "count") + 
	scale_fill_brewer(palette = "RdPu") + 
	theme_bw() +
	labs(y = "N", fill = "Self-Assessed Likelihood of\nJob Being Automated") +
	theme(axis.text.x = element_text(angle=45, hjust = 1))
dev.off()

####################################################



####################################################
# Figure 4. Partisanship and vulnerability
data$Party <- rep(NA, length(data$PID_BIN))
data$Party[data$PID_BIN == -1] <- "Republican"
data$Party[data$PID_BIN == 0] <- "Independent"
data$Party[data$PID_BIN == 1] <- "Democrat"

data$YourJobFear <- rep(NA, length(data$PID_BIN))
data$YourJobFear[data$YOUR_JOB_AUTO == 1] <- "Not At All Likely"
data$YourJobFear[data$YOUR_JOB_AUTO == 2] <- "Not Very Likely"
data$YourJobFear[data$YOUR_JOB_AUTO == 3] <- "Somewhat Likely"
data$YourJobFear[data$YOUR_JOB_AUTO == 4] <- "Very Likely"

fear.plot <- na.omit(data.frame(Vulnerability = data$ECON_VULNERABLE,
				Party = data$Party, YourJobFear = data$YourJobFear))
				
# Figure 4. Who feels vulnerable?
ggplot(data = fear.plot, aes(x = Vulnerability, fill = YourJobFear, lty = YourJobFear)) +
	geom_density(alpha = 0.4) +
	facet_grid(rows = vars(Party)) +
	scale_fill_brewer(palette = "RdPu") + 
	labs(y = "Density", x = "Self-Rated Economic Vulnerability\n(0 = Least Vulnerable, 100 = Most Vulnerable)", 
		lty = "Self-Assessed Likelihood of\nOwn Job Being Automated", fill = "Self-Assessed Likelihood of\nOwn Job Being Automated") +
	theme_bw()
dev.off()

####################################################







####################################################
# Table 1. Simple, based on just threat of automation.
#  Column 1. All respondents
model.simpleindus.threat <- lm(YOUR_JOB_AUTO ~ as.factor(YOUR_JOB_AUTO_LOW_MID_HIGH), data = data)
summary(model.simpleindus.threat)

#  Column 2. Just Republicans
model.simpleindus.threat.reps <- lm(YOUR_JOB_AUTO ~ as.factor(YOUR_JOB_AUTO_LOW_MID_HIGH), 
	data = data, subset = PID_BIN < 0)
summary(model.simpleindus.threat.reps)

#  Column 3. Just Democrats
model.simpleindus.threat.dems <- lm(YOUR_JOB_AUTO ~ as.factor(YOUR_JOB_AUTO_LOW_MID_HIGH), 
	data = data, subset = PID_BIN > 0)
summary(model.simpleindus.threat.dems)

model.simpleindus.coef <- c("Intercept", "Medium Threat", "High Threat", "Freelance", 
	"Self-Employed")

model.simpleindus.table <- NULL # initialize

for(i in 1:length(model.simpleindus.coef)) {
	model.simpleindus.table <- rbind(model.simpleindus.table, 
	paste0(model.simpleindus.coef[i], " & ", round(coef(model.simpleindus.threat)[i], digits = 3),
		" & ", round(coef(model.simpleindus.threat.reps)[i], digits = 3),
		" & ", round(coef(model.simpleindus.threat.dems)[i], digits = 3), " \\"),
	paste0( " ", " & (", round(sqrt(diag(vcov(model.simpleindus.threat)))[i], digits = 3),
		") & (", round(sqrt(diag(vcov(model.simpleindus.threat.reps)))[i], digits = 3),
		") & (", round(sqrt(diag(vcov(model.simpleindus.threat.dems)))[i], digits = 3), ") \\"))
}		

model.simpleindus.table <- rbind(model.simpleindus.table, 
		paste0("$R^2$", " & ", round(summary(model.simpleindus.threat)$adj.r.square, digits = 3),
		 " & ", round(summary(model.simpleindus.threat.reps)$adj.r.square, digits = 3),
		 " & ", round(summary(model.simpleindus.threat.dems)$adj.r.square, digits = 3), " \\"),
		paste0("$N$", " & ", length(residuals(model.simpleindus.threat)),
		 " & ", length(residuals(model.simpleindus.threat.reps)),
		 " & ", length(residuals(model.simpleindus.threat.dems)), " \\"))

model.simpleindus.table

####################################################





####################################################
# Table 2. Simple, based on just threat of automation x fear of automation
#  Column 1. All respondents
model.ubi.indus.threat <- lm(UBI_SUPPORT ~ as.factor(YOUR_JOB_AUTO_LOW_MID_HIGH)*YOUR_JOB_AUTO, 
	data = data)
summary(model.ubi.indus.threat)

#  Column 2. Just Republicans
model.ubi.indus.threat.rep <- lm(UBI_SUPPORT ~ as.factor(YOUR_JOB_AUTO_LOW_MID_HIGH)*YOUR_JOB_AUTO, 
	data = data, subset = PID_BIN < 0)
summary(model.ubi.indus.threat.rep)

#  Column 3. Just Democrats
model.ubi.indus.threat.dem <- lm(UBI_SUPPORT ~ as.factor(YOUR_JOB_AUTO_LOW_MID_HIGH)*YOUR_JOB_AUTO, 
	data = data, subset = PID_BIN > 0)
summary(model.ubi.indus.threat.dem)


model.ubi.indus.coef <- c("Intercept", "Medium Threat", "High Threat", "Freelance", 
	"Self-Employed", "Fear (of Your Job Automated)", "Medium:Fear", "High:Fear", 
	"Freelance:Fear", "Self:Fear")

model.ubi.indus.table <- NULL # initialize

for(i in 1:length(model.ubi.indus.coef)) {
	model.ubi.indus.table <- rbind(model.ubi.indus.table, 
	paste0(model.ubi.indus.coef[i], " & ", round(coef(model.ubi.indus.threat)[i], digits = 3),
		" & ", round(coef(model.ubi.indus.threat.rep)[i], digits = 3),
		" & ", round(coef(model.ubi.indus.threat.dem)[i], digits = 3), " \\"),
	paste0( " ", " & (", round(sqrt(diag(vcov(model.ubi.indus.threat)))[i], digits = 3),
		") & (", round(sqrt(diag(vcov(model.ubi.indus.threat.rep)))[i], digits = 3),
		") & (", round(sqrt(diag(vcov(model.ubi.indus.threat.dem)))[i], digits = 3), ") \\"))
}		

model.ubi.indus.table <- rbind(model.ubi.indus.table, 
		paste0("$R^2$", " & ", round(summary(model.ubi.indus.threat)$adj.r.square, digits = 3),
		 " & ", round(summary(model.ubi.indus.threat.rep)$adj.r.square, digits = 3),
		 " & ", round(summary(model.ubi.indus.threat.dem)$adj.r.square, digits = 3), " \\"),
		paste0("$N$", " & ", length(residuals(model.ubi.indus.threat)),
		 " & ", length(residuals(model.ubi.indus.threat.rep)),
		 " & ", length(residuals(model.ubi.indus.threat.dem)), " \\"))

model.ubi.indus.table

####################################################







####################################################
# Figure 5. Fit the predictions from Table 2.

prediction.data <- expand.grid(YOUR_JOB_AUTO = 1:4, 
							YOUR_JOB_AUTO_LOW_MID_HIGH = levels(as.factor(data$YOUR_JOB_AUTO_LOW_MID_HIGH)))

rep.predictions <- data.frame(Fear = prediction.data$YOUR_JOB_AUTO,
							UBI.Support = predict(model.ubi.indus.threat.rep, prediction.data),
							Party = "Republican",
							threatnum = prediction.data$YOUR_JOB_AUTO_LOW_MID_HIGH)

dem.predictions <- data.frame(Fear = prediction.data$YOUR_JOB_AUTO,
							UBI.Support = predict(model.ubi.indus.threat.dem, prediction.data),
							Party = "Democrat",
							threatnum = prediction.data$YOUR_JOB_AUTO_LOW_MID_HIGH)

ubi.gg <- rbind(rep.predictions, dem.predictions)

ubi.gg$Threat <- recode(ubi.gg$threatnum, `1` = "Low Threat", `2` = "Medium Threat", `3` = "High Threat",
											`4` = "Freelance", `5` = "Self-Employed")



# Figure 5. Support by threat and party (Figure 5)
ggplot(data = ubi.gg, aes(x = Fear, y = UBI.Support, col = Threat)) +
	geom_line(lwd = 2) + 
	scale_color_brewer(palette = "RdPu") + 
	facet_wrap(~Party) + 
	labs(x = "Self-Assessed Likelihood of Own Job Being Automated\n(1 = Not at all likely, 4 = Very likely)", 
		y = "Support for UBI (1 = Not At All Likely to Support,\n5 = Very Likely to Support)",
		col = "Actual Threat\nof Automation") +
	theme_bw()
dev.off()

####################################################

























####################################################
# Appendix Table 1. Prediction of UBI with controls
#  Column 1. All respondents
model.ubi.indus.threat.controls.noint <- lm(UBI_SUPPORT ~
	as.factor(YOUR_JOB_AUTO_LOW_MID_HIGH) + YOUR_JOB_AUTO + 
	DEFICIT + INCOME + EDUC + as.factor(TREATMENT) + IDEOLOGY + FEMALE + RETRO + PROSPECT + 
	ECON_VULNERABLE + VOTEYES, 
	data = data)
summary(model.ubi.indus.threat.controls.noint)

#  Column 2. Just Republicans
model.ubi.indus.threat.controls.noint.dem <- lm(UBI_SUPPORT ~
	as.factor(YOUR_JOB_AUTO_LOW_MID_HIGH) + YOUR_JOB_AUTO + 
	DEFICIT + INCOME + EDUC + as.factor(TREATMENT) + IDEOLOGY + FEMALE + RETRO + PROSPECT + 
	ECON_VULNERABLE + VOTEYES, 
	data = data, subset = PID_BIN > 0)
summary(model.ubi.indus.threat.controls.noint.dem)

#  Column 3. Just Democrats
model.ubi.indus.threat.controls.noint.rep <- lm(UBI_SUPPORT ~
	as.factor(YOUR_JOB_AUTO_LOW_MID_HIGH) + YOUR_JOB_AUTO + 
	DEFICIT + INCOME + EDUC + as.factor(TREATMENT) + IDEOLOGY + FEMALE + RETRO + PROSPECT + 
	ECON_VULNERABLE + VOTEYES, 
	data = data, subset = PID_BIN < 0)
summary(model.ubi.indus.threat.controls.noint.rep)


model.ubi.indus.controls.coef <- c("Intercept", "Medium Threat", "High Threat", "Freelance", 
	"Self-Employed", "Fear (of Your Job Automated)", "Concern for Deficit", "Income", 
	"Education", "UBIF", "UBIF", "UBIF", "UBIF", "UBIF", "UBIF", "UBIF", "UBIF", 
	"Ideology", "Female", "Retrospections", "Prospections", "Economic Vulnerability", 
	"Voter")

model.ubi.indus.controls.table <- NULL # initialize

for(i in 1:length(model.ubi.indus.controls.coef)) {
	model.ubi.indus.controls.table <- rbind(model.ubi.indus.controls.table, 
	paste0(model.ubi.indus.controls.coef[i], " & ", 
		round(coef(model.ubi.indus.threat.controls.noint)[i], digits = 3),
		" & ", round(coef(model.ubi.indus.threat.controls.noint.rep)[i], digits = 3),
		" & ", round(coef(model.ubi.indus.threat.controls.noint.dem)[i], digits = 3), " \\"),
	paste0( " ", " & (", round(sqrt(diag(vcov(model.ubi.indus.threat.controls.noint)))[i], digits = 3),
		") & (", round(sqrt(diag(vcov(model.ubi.indus.threat.controls.noint.rep)))[i], digits = 3),
		") & (", round(sqrt(diag(vcov(model.ubi.indus.threat.controls.noint.dem)))[i], digits = 3), ") \\"))
}		

model.ubi.indus.controls.table <- rbind(model.ubi.indus.controls.table, 
		paste0("$R^2$", " & ", round(summary(model.ubi.indus.threat.controls.noint)$adj.r.square, digits = 3),
		 " & ", round(summary(model.ubi.indus.threat.controls.noint.rep)$adj.r.square, digits = 3),
		 " & ", round(summary(model.ubi.indus.threat.controls.noint.dem)$adj.r.square, digits = 3), " \\"),
		paste0("$N$", " & ", length(residuals(model.ubi.indus.threat.controls.noint)),
		 " & ", length(residuals(model.ubi.indus.threat.controls.noint.rep)),
		 " & ", length(residuals(model.ubi.indus.threat.controls.noint.dem)), " \\"))

model.ubi.indus.controls.table
####################################################














